***** Programmer: Adaobi Anakwe
***** Analysis: Trend Analysis in Men's Preconception Health: Variation by race]
***** Start date: August 23, 2021
****************************************************************

*****************************************************************
****Importing and cleaning data set (labeling variables)******
*****************************************************************

log using "D:\MYFILES\NSFG\data_management\Trends1", replace 
import sas using "D:\MYFILES\NSFG\combine1119.sas7bdat"
describe
summarize
*********************************************************Examining missingness on all variables **********
*****************************************************
tab gender, m 
tab pir, m 
tab educ, m 
tab pregint, m  
tab insurance, m 
tab reprohlth, m
tab racial, m 
tab immigration, m 
tab marital, m 
tab sexprt2, m 
tab sexrisk2, m 
tab condomuse, m 
tab genhealth, m 
tab alcuse, m 
tab druguse, m 
tab sti2, m 
tab bmicat2, m  
tab msa, m 
tab workcons2, m 
tab SMK100, m 
tab SMOKE12, m
tab PARITY, m
tab year, m
tab HADSEX, m
tab VRY1STAG, m

****************************************************
/************Coding and labeling variables******/
****************************************************
label var sexprt2 "Number of sexual partners"
label define sexprt2 0 "none to 1" 1 "more than 1"
label values sexprt2 sexprt2
tab sexprt2, m

label var sexrisk2 "Sexual risk behavior"
label define sexrisk2 0 "No risk" 1 "high risk"
label values sexrisk2 sexrisk2
tab sexrisk2, m

label var condomuse "Condom use"
label define condomuse 0 "consistent use" 1 "inconsistent use to no use"
label values  condomuse condomuse 
tab condomuse, m 

label var genhealth "General health"
label define genhealth 0 "'excellent to good" 1 "fair to poor"
label values genhealth	genhealth
tab genhealth, m 

label var alcuse "Alcohol use"
label define alcuse 0 "no to low alcohol use" 1 "medium to high alcohol use"
label values alcuse alcuse
tab alcuse, m

label var druguse "Past 12 months drug use"
label define druguse 0 "no" 1 "yes"
label values druguse druguse
tab druguse, m

label var sti2 "Sexual Transmitted diseases"
label define sti2 0 "no" 1 "yes"
label values sti2 sti2
tab sti2, m
 
label var bmicat2 "Body Mass Index"
label define bmicat2 0 "underweight to normal weight" 1 "overweight to obese"
label value bmicat2 bmicat2
tab bmicat2
tab bmicat2, m 

label var pregint "Pregnancy intentions"
label define pregint 0 "no intention" 1 "intends pregnancy" 2 "ambivalent"
label values pregint pregint
tab pregint, m

label var racial "race/ethnicity"
label define racial 0  "nh white" 1 "nh black" 2 "hispanic" 3 "other"
label values racial racial
tab racial, m

label var pir "poverty income ratio"
label define pir 0 "below 100%" 1 "at 100% and above"
label values pir pir
tab pir, m

label var educ "educational attainment"
label define educ 0 "less than HS" 1 "high school" 2 "some college" 3 "college or more"
label values educ educ
tab educ, m

label var insurance "insurance status in the last 12 months"
label define insurance 0 "no insurance" 1 "had insurance"
label values insurance insurance
tab insurance, m
 
label var reprohlth "Reproductive health visit"
label define reprohlth 0 "No visit" 1 "non reproductive visit" 2 "reproductive health visit"
label values reprohlth reprohlth
tab reprohlth, m

label var immigration "immigration status"
label define immigration 0 "non-immigrant" 1 "immigrant"
label values immigration immigration
tab immigration, m

label var marital "marital status"
label define marital 0 "married" 1 "separated/widowed/divorced" 2 "never married"
label values marital marital
tab marital, m

label var msa "metropolitan area"
label define msa 0 "urban" 1 "suburban" 2 "rural"
label values msa msa
tab msa, m

label var workcons2 "employment consistency"
label define workcons2 0 "not employed" 1 "unstable employment" 2 "stable employment"
label values workcons2 workcons2
tab workcons2, m

label var gender "male vs female"
label define gender 0 "male" 1 "female"
label values gender gender
tab gender, m

recode AGER 20/24=0 25/29=1 30/34=2 35/39=3 40/44=4, generate (agecat)
gen agecat2=.
replace agecat2 = 0 if agecat == 0
replace agecat2 = 1 if agecat == 1
replace agecat2 = 2 if agecat == 2
replace agecat2 = 3 if agecat == 3
replace agecat2 = 4 if agecat == 4
label var agecat2 "Participants age categories"
label define agecat2 0 "20 to 24 years" 1 "25 to 29 years" 2 "30 to 34 years" 3 "35 to 39 years" 4 "40 to 44 years"
label values agecat2 agecat2
tab agecat2, m

tab year

/*** generating the year variable******/
gen yeara = . 
replace yeara = 0 if year == 2012 
replace yeara = 1 if year == 2014
replace yeara = 2 if year == 2016
replace yeara = 3 if year == 2018
label var yeara "Midpoint Year"
label define yeara 0 "2012" 1 "2014" 2 "2016" 3 "2018"
label values yeara yeara 
tab yeara 

gen year2 = yeara*yeara
label var year2 "Year, Quadratic"
gen year3 = yeara*yeara*yeara
label var year3 "Year, Cubic"

tab year2
tab year3
tab gender

generate pchscore= sexprt2+sexrisk2+condomuse+alcuse+druguse+sti2+genhealth+bmicat2
label var pchscore "preconception health score"
tab pchscore, m

recode pchscore 0=0 1=1 2=2 3=3 4=4 5=4 6=4 7=4 8=4, generate (pchord)
gen pchord2 =.
replace pchord2= 0 if pchord==0
replace pchord2=1 if pchord==1
replace pchord2=2 if pchord==2
replace pchord2=3 if pchord==3
replace pchord2=4 if pchord==4
label var pchord2 "preconception risk ordered"
label define pchord2 0 "no risk factor" 1 "one risk factor" 2 " two risk factors" 3 "three risk factors" 4 " 4 or more risk factors"
label values pchord2 pchord2
tab pchord2, m

tab pchord2 pchscore
tab pchord2 if gender==0, m

/***********Missingness check**********************************************/
***************************************************************************
misschk pchscore gender sexprt2 sexrisk2 condomuse alcuse druguse sti2  genhealth bmicat2 pregint racial pir educ insurance immigration marital msa workcons2 AGER BIOKIDS VRY1STAG year if gender !=., gen(cov_miss_)

misschk pchscore gender sexprt2 sexrisk2 condomuse alcuse druguse sti2  genhealth bmicat2 pregint racial pir educ insurance reprohlth immigration marital msa workcons2 AGER BIOKIDS VRY1STAG year if gender == 0, gen(miss)

*****************************************************************************************************
******** 6.1% missingness on sexrisk2 and 6.2% missingness on condomuse******************************
******** A missingness variable will be created and checked for significant association*************
******** Since missingness less than 10%, will use listwise deletion in regression models***********
******************************************************************************************
log close


log using "D:\MYFILES\NSFG\data_management\descriptives", replace 

/************* Declaring the survey weight******/

svyset panelvar [pweight=WEIGHTVAR], strata(stratvar) vce(linearized) singleunit(centered)

global studypop  "if gender==0" /***All races***/
global studypop1 "if gender==0 & racial==0" /***NH White men***/
global studypop2 "if gender==0 & racial==1" /***NH Black men***/
global studypop3 "if gender==0 & racial==2" /*** Hispanic men***/


/*********************RQ1: Is there a linear, cubic or quadriatric trend in men's preconception health risk factors******************************/
******************************************************************************************
******************************************************************************************
******************************************************************************************

******* Descriptives overall population [take output for only totals]********************

foreach x in sexprt2 sexrisk2 condomuse alcuse druguse sti2  genhealth bmicat2 pchscore pregint racial pir educ insurance reprohlth immigration marital msa workcons2 agecat2 {
	estpost svy, subpop($studypop): tab `x' year, col percent se 
	estimates store  tab_`x' 
				} 
svy linearized, subpop(if gender==0) : mean BIOKIDS AGER VRY1STAG, over(year)
	estat sd
svy linearized, subpop(if gender==0) : mean BIOKIDS AGER VRY1STAG pchscore
	estat sd
**********Descriptive for NH White men*********
foreach x in sexprt2 sexrisk2 condomuse alcuse druguse sti2  genhealth bmicat2 pchscore pregint  pir educ insurance reprohlth immigration marital msa workcons2 agecat2 {
	estpost svy, subpop($studypop1): tab `x' year, col percent se 
	estimates store  tab_`x' 
				} 
svy linearized, subpop(if gender==0 & racial==0) : mean BIOKIDS AGER VRY1STAG, over(year)
	estat sd
svy linearized, subpop(if gender==0 & racial==0) : mean BIOKIDS AGER VRY1STAG pchscore
	estat sd

******* Descriptives for NH Black men************
foreach x in sexprt2 sexrisk2 condomuse alcuse druguse sti2 genhealth bmicat2 pchscore pregint  pir educ insurance reprohlth immigration marital msa workcons2 agecat2 {
	estpost svy, subpop($studypop2): tab `x' year, col percent se 
	estimates store  tab_`x' 
				} 
svy linearized, subpop(if gender==0 & racial==1) : mean BIOKIDS AGER VRY1STAG, over(year)
	estat sd
svy linearized, subpop(if gender==0 & racial==1) : mean BIOKIDS AGER VRY1STAG pchscore
	estat sd
****** Descriptives for Hispanic**************

foreach x in sexprt2 sexrisk2 condomuse alcuse druguse sti2  genhealth bmicat2 pchscore pregint  pir educ insurance reprohlth immigration marital msa workcons2 agecat2 {
	estpost svy, subpop($studypop3): tab `x' year, col percent se 
	estimates store  tab_`x' 
				} 
svy linearized, subpop(if gender==0 & racial==2) : mean BIOKIDS AGER VRY1STAG, over(year)
	estat sd
svy linearized, subpop(if gender==0 & racial==2) : mean BIOKIDS AGER VRY1STAG pchscore
	estat sd

log close

 ***************Exporting results to .csv***************
*foreach x in sexprt2 sexrisk2 condomuse alcuse druguse sti2 genhealth bmicat2 pregint racial pir educ insurance immigration marital msa workcons2 agecat2 { 
*	esttab tab_`x' using "C:\Users\Owner\Desktop\MYFILES\NSFG\descriptives.csv", nomtitles noisily cells(cell(fmt(3))lb ub) label ///
*		varlabels (`e(labels)') eqlabels(`e(eqlabels)') ///
*		nonum unstack onecell noobs baselevels ///
*		append stats(p_Pear N_sub, fmt(3 0) label ("P-value")) drop(Total)
*}

log using "D:\MYFILES\NSFG\data_management\TrendsTest_overall", replace

****************************************************************************************
**************************TEST FOR TRENDS****************************************
****************************************************************************************

svy linearized, subpop(if gender==0): logit sexprt2 i.yeara, eform
contrast p.yeara, noeffects 
svy linearized, subpop(if gender==0): logit sexrisk2 i.yeara, eform
contrast p.yeara, noeffects
svy linearized, subpop(if gender==0): logit condomuse i.yeara, eform
contrast p.yeara, noeffects
svy linearized, subpop(if gender==0): logit alcuse i.yeara, eform
contrast p.yeara, noeffects
svy linearized, subpop(if gender==0): logit druguse i.yeara, eform
contrast p.yeara, noeffects
svy linearized, subpop(if gender==0): logit sti2 i.yeara, eform
contrast p.yeara, noeffects
svy linearized, subpop(if gender==0): logit genhealth i.yeara, eform
contrast p.yeara, noeffects
svy linearized, subpop(if gender==0): logit bmicat2 i.yeara, eform
contrast p.yeara, noeffects
svy linearized, subpop(if gender==0): logit pir i.yeara, eform
contrast p.yeara, noeffects
svy linearized, subpop(if gender==0): logit insurance i.yeara, eform
contrast p.yeara, noeffects
svy linearized, subpop(if gender==0): logit reprohlth i.yeara, eform
contrast p.yeara, noeffects
svy linearized, subpop(if gender==0): logit educ i.yeara, eform
contrast p.yeara, noeffects
svy linearized, subpop(if gender==0): logit immigration i.yeara, eform
contrast p.yeara, noeffects
svy linearized, subpop(if gender==0): logit marital i.yeara, eform
contrast p.yeara, noeffects
svy linearized, subpop(if gender==0): logit msa i.yeara, eform
contrast p.yeara, noeffects
svy linearized, subpop(if gender==0): logit workcons2 i.yeara, eform
contrast p.yeara, noeffects
svy linearized, subpop(if gender==0): logit agecat2 i.yeara, eform
contrast p.yeara, noeffects
svy linearized, subpop(if gender==0): regress pchscore i.yeara 
****************************************************************************************
*********************EXAMINING TRENDS BY RACE/ETHNICITY*********************************
****************************************************************************************
svy linearized, subpop(if racial==0 & gender==0): logit sexprt2 i.yeara, eform
contrast p.yeara, noeffects 
svy linearized, subpop(if racial==1 & gender==0): logit sexprt2 i.yeara, eform
contrast p.yeara, noeffects 
svy linearized, subpop(if racial==2 & gender==0): logit sexprt2 i.yeara, eform
contrast p.yeara, noeffects 
 
svy linearized, subpop(if racial==0 & gender==0): logit sexrisk2 i.yeara, eform
contrast p.yeara, noeffect
svy linearized, subpop(if racial==1 & gender==0): logit sexrisk2 i.yeara, eform
contrast p.yeara, noeffect
svy linearized, subpop(if racial==2 & gender==0): logit sexrisk2 i.yeara, eform
contrast p.yeara, noeffect

svy linearized, subpop(if racial==0 & gender==0): logit condomuse i.yeara, eform
contrast p.yeara, noeffect
svy linearized, subpop(if racial==1 & gender==0): logit condomuse i.yeara, eform
contrast p.yeara, noeffect
svy linearized, subpop(if racial==2 & gender==0): logit condomuse i.yeara, eform
contrast p.yeara, noeffect

svy linearized, subpop(if racial==0 & gender==0): logit alcuse i.yeara, eform
contrast p.yeara, noeffect
svy linearized, subpop(if racial==1 & gender==0): logit alcuse i.yeara, eform
contrast p.yeara, noeffect
svy linearized, subpop(if racial==2 & gender==0): logit alcuse i.yeara, eform
contrast p.yeara, noeffect

svy linearized, subpop(if racial==0 & gender==0): logit druguse i.yeara, eform
contrast p.yeara, noeffect
svy linearized, subpop(if racial==1 & gender==0): logit druguse i.yeara, eform
contrast p.yeara, noeffect
svy linearized, subpop(if racial==2 & gender==0): logit druguse i.yeara, eform
contrast p.yeara, noeffect

svy linearized, subpop(if racial==0 & gender==0): logit sti2 i.yeara, eform
contrast p.yeara, noeffect
svy linearized, subpop(if racial==1 & gender==0): logit sti2 i.yeara, eform
contrast p.yeara, noeffect
svy linearized, subpop(if racial==2 & gender==0): logit sti2 i.yeara, eform
contrast p.yeara, noeffect

svy linearized, subpop(if racial==0 & gender==0): logit genhealth i.yeara, eform
contrast p.yeara, noeffect
svy linearized, subpop(if racial==1 & gender==0): logit genhealth i.yeara, eform
contrast p.yeara, noeffect
svy linearized, subpop(if racial==2 & gender==0): logit genhealth i.yeara, eform
contrast p.yeara, noeffect

svy linearized, subpop(if racial==0 & gender==0): logit bmicat2 i.yeara, eform
contrast p.yeara, noeffect
svy linearized, subpop(if racial==1 & gender==0): logit bmicat2 i.yeara, eform
contrast p.yeara, noeffect
svy linearized, subpop(if racial==2 & gender==0): logit bmicat2 i.yeara, eform
contrast p.yeara, noeffect

svy linearized, subpop(if racial==0 & gender==0): regress pchscore i.yeara 
svy linearized, subpop(if racial==1 & gender==0): regress pchscore i.yeara 
svy linearized, subpop(if racial==2 & gender==0): regress pchscore i.yeara 

********************************************************************************************
******** STEP 3A: Examining associations between pchscore & racial and demographic variables
********************************************************************************************
svy linearized, subpop(if gender==0): regress pchscore i.pir
svy linearized, subpop(if gender==0): regress pchscore i.insurance
svy linearized, subpop(if gender==0): regress pchscore i.educ
svy linearized, subpop(if gender==0): regress pchscore i.immigration
svy linearized, subpop(if gender==0): regress pchscore i.marital
svy linearized, subpop(if gender==0): regress pchscore i.msa
svy linearized, subpop(if gender==0): regress pchscore i.workcons2
svy linearized, subpop(if gender==0): regress pchscore i.agecat2
svy linearized, subpop(if gender==0): regress pchscore i.reprohlth
svy linearized, subpop(if gender==0): regress pchscore BIOKIDS
svy linearized, subpop(if gender==0): regress pchscore VRY1STAG

svy linearized, subpop(if gender==0): mlogit racial i.pir, rrr
svy linearized, subpop(if gender==0): mlogit racial i.insurance, rrr
svy linearized, subpop(if gender==0): mlogit racial i.educ, rrr
svy linearized, subpop(if gender==0): mlogit racial i.immigration, rrr
svy linearized, subpop(if gender==0): mlogit racial i.marital, rrr
svy linearized, subpop(if gender==0): mlogit racial i.msa, rrr
svy linearized, subpop(if gender==0): mlogit racial i.workcons2, rrr
svy linearized, subpop(if gender==0): mlogit racial i.agecat2, rrr
svy linearized, subpop(if gender==0): mlogit racial BIOKIDS
svy linearized, subpop(if gender==0): mlogit racial VRY1STAG
svy linearized, subpop(if gender==0): mlogit racial i.reprohlth, rrr

*******************STEP 3B:Checking the linear relationship assumption between exposure (racial) and outcome (pchscore)********/

svy linearized, subpop(if gender==0): regress pchscore i.racial

rvfplot, recast(scatter)
rvfplot, yline(0)

twoway (scatter pchscore racial if gender==0)
****************************************************************************************
********************STEP 3C: Regression**************************************************
****************************************************************************************
svy linearized, subpop(if gender==0): regress pchscore i.racial 
svy linearized, subpop(if gender==0): regress pchscore i.yeara 
svy linearized, subpop(if gender==0): regress pchscore i.racial i.year i.racial##i.yeara
testparm i.racial##i.yeara /****since statistically significant interaction term, results will be strattified by race*******/

**************************
**multivariate linear regression- full sample
**************************
svy linearized, subpop(if gender==0): regress pchscore i.racial i.year i.pir i.insurance i.educ i.immigration i.marital i.msa i.workcons2 i.agecat2 i.reprohlth BIOKIDS VRY1STAG

svy linearized, subpop(if gender==0): regress pchscore i.racial i.year i.pir i.insurance i.immigration i.marital i.msa i.agecat2 i.reprohlth VRY1STAG

****************************************
*multivariate linear regression- subsample: NH WHITE
*****************************************
svy linearized, subpop(if racial==0 & gender==0): regress pchscore i.year i.pir i.insurance i.educ i.immigration i.marital i.msa i.workcons2 i.agecat2 i.reprohlth BIOKIDS VRY1STAG

svy linearized, subpop(if racial==0 & gender==0): regress pchscore i.year i.insurance i.marital i.msa i.reprohlth VRY1STAG

****************************************
*multivariate linear regression- subsample: NH BLACK
*****************************************
svy linearized, subpop(if racial==1 & gender==0): regress pchscore i.year i.pir i.insurance i.educ i.immigration i.marital i.msa i.workcons2 i.agecat2 i.reprohlth BIOKIDS VRY1STAG

svy linearized, subpop(if racial==1 & gender==0): regress pchscore i.year i.marital i.reprohlth  VRY1STAG

****************************************
*multivariate linear regression- subsample: Hispanic
*****************************************
svy linearized, subpop(if racial==2 & gender==0): regress pchscore i.year i.pir i.insurance i.educ i.immigration i.marital i.msa i.workcons2 i.agecat2 i.reprohlth BIOKIDS VRY1STAG

svy linearized, subpop(if racial==2 & gender==0): regress pchscore i.year i.pir i.immigration i.marital i.agecat2 i.reprohlth VRY1STAG

************************************************************************************ GENERALIZED ORDINARY LOGIT MODEL***************************
*****************************************************************************

********STEP 1: Constructing the order preconception variable**************

recode pchscore 0=0 1=1 2=2 3=3 4=4 5=4 6=4 7=4 8=4, generate (pchord)
gen pchord2 =.
replace pchord2= 0 if pchord==0
replace pchord2=1 if pchord==1
replace pchord2=2 if pchord==2
replace pchord2=3 if pchord==3
replace pchord2=4 if pchord==4
label var pchord2 "preconception risk ordered"
label define pchord2 0 "no risk factor" 1 "one risk factor" 2 " two risk factors" 3 "three risk factors" 4 " 4 or more risk factors"
label values pchord2 pchord2
tab pchord2, m

tab pchord2 pchscore
tab pchord2 if gender==0, m

********* STEP 1b: Testing the proportional odds assumption

svy linearized, subpop(if gender==0): gologit2 pchord2 i.racial i.year, pl

svy linearized, subpop(if gender==0): gologit2 pchord2 i.racial i.year, npl

svy linearized, subpop(if gender==0): ologit pchord2 i.racial i.year  

svy linearized, subpop(if gender==0): gologit2 pchord2 i.racial i.year, auto(.01) lrf

**********************

gsvy: gologit2 pchord2 i.racial i.yeara i.pir i.insurance i.educ i.immigration i.marital i.msa i.workcons2 i.agecat2 i.reprohlth BIOKIDS VRY1STAG, auto(.025)

gsvy: gologit2 pchord2 i.racial i.yeara i.pir i.insurance i.educ i.immigration i.marital i.msa i.workcons2 i.agecat2 i.reprohlth BIOKIDS VRY1STAG, auto(.025)
gsvy: gologit2 pchord2 i.racial i.yeara i.pir i.insurance i.educ i.immigration i.marital i.msa i.workcons2 i.agecat2 i.reprohlth BIOKIDS VRY1STAG, pl(0b.reprohlth 1.reprohlth 2.reprohlth 0b.racial 1.racial 2.racial 3.racial 0b.yeara 1.yeara 2.yeara 3.yeara 0b.pir 1.pir 0b.insurance 1.insurance 0b.educ 1.educ 2.educ 0b.immigration 0b.marital 0b.msa 1.msa 2.msa 0b.workcons2 1.workcons2 2.workcons2 0b.agecat2 VRY1STAG)
mchange, amount(margin)


***************************************************************
*log close
**********
**********To calculate the means with p-values (code from Kasim)*******
/*	Continuous Covariates */
foreach x in oschlyrs r8agey_m r8bmi r8conden{	
	*Run means and SD to create standardized scores
	svy, subpop($study_subp): mean `x'
	estat sd
	return list 
	mat `x'_mean_sd_1 = (r(mean)\r(sd))
	*Run means and SD over Diabetes Status for standardized scores
	svy, subpop($study_subp): mean `x', over(diab3k)
	estat sd
	return list
	mat `x'_mean_sd_2 = (r(mean)\r(sd))
	*Creates a combined matrix
	mat `x'_mean_sd = (`x'_mean_sd_2, `x'_mean_sd_1)
	mat list `x'_mean_sd
	mat rownames `x'_mean_sd = "Mean_`x'" "SD_`x'"
	svy, subpop($study_subp): reg `x' i.diab3k
	estimate store out_`x'
	mat `x'_p_value = e(p)
	mat `x'_subpop = e(N_sub)
	mat rownames `x'_p_value = P-Value
	mat colnames `x'_p_value = ""
	mat rownames `x'_subpop = N_Sub
	mat colnames `x'_subpop = ""
	}

*********** To export in as CSV in Excel****************************
	/*	Continuous Covariates */

foreach x in oschlyrs r8agey_m r8bmi r8conden{	
	esttab matrix(`x'_mean_sd, fmt(2)) using descriptives.csv, nomtitles ///
	label append nonum ///
	collabels("Low" "Moderate" "High")
	esttab matrix(`x'_p_value, fmt(2)) using descriptives.csv, nomtitles collabels(none) ///
	label append nonum
	esttab matrix(`x'_subpop, fmt(2)) using descriptives.csv, nomtitles collabels(none) ///
	label append nonum
}
